Instructions:

• Load the boundaries of Danish municipalities from data/ folder, convert to simple feature and transform to CRS 25832.
• Sort the NAME_2 field to see how the Danish municipalities are spelled. You may need to change them later for the spatial data to join the attributes.

# Load the spatial data, project to UTM
mun_sp<- readRDS(________) # it's the gadm_... 2.rds dataset
mun_sf <- ________(mun_sp)
mun <- st_transform(_____, crs = ____)

# Plot so as to check correct location and complete coverage


# Check the names


# Straighten the names (return here after Task 2)

Simple feature collection with 99 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 441745.6 ymin: 6049775 xmax: 892801.1 ymax: 6402207
Projected CRS: ETRS89 / UTM zone 32N
First 10 features:
      GID_0  NAME_0   GID_1      NAME_1 NL_NAME_1     GID_2      NAME_2
36829   DNK Denmark DNK.1_1 Hovedstaden      <NA> DNK.1.1_1 Albertslund
36665   DNK Denmark DNK.1_1 Hovedstaden      <NA> DNK.1.2_1     Allerød
36779   DNK Denmark DNK.1_1 Hovedstaden      <NA> DNK.1.3_1    Ballerup
37011   DNK Denmark DNK.1_1 Hovedstaden      <NA> DNK.1.4_1    Bornholm
36901   DNK Denmark DNK.1_1 Hovedstaden      <NA> DNK.1.5_1     Brøndby
      VARNAME_2 NL_NAME_2  TYPE_2    ENGTYPE_2 CC_2   HASC_2
36829      <NA>      <NA> Kommune Municipality <NA> DK.HS.AB
36665      <NA>      <NA> Kommune Municipality <NA> DK.HS.AL
36779      <NA>      <NA> Kommune Municipality <NA> DK.HS.BA
37011      <NA>      <NA> Kommune Municipality <NA> DK.HS.BO
36901      <NA>      <NA> Kommune Municipality <NA> DK.HS.BR
                            geometry
36829 MULTIPOLYGON (((712057 6173...
36665 MULTIPOLYGON (((700891 6191...
36779 MULTIPOLYGON (((715156 6178...
37011 MULTIPOLYGON (((878103.7 61...
36901 MULTIPOLYGON (((716929 6168...
 [ reached 'max' / getOption("max.print") -- omitted 5 rows ]

 [1] "Albertslund"       "Allerød"           "Ballerup"         
 [4] "Bornholm"          "Brøndby"           "Christiansø"      
 [7] "Dragør"            "Egedal"            "Fredensborg"      
[10] "Frederiksberg"     "Frederikssund"     "Furesø"           
[13] "Gentofte"          "Gladsaxe"          "Glostrup"         
[16] "Gribskov"          "Halsnæs"           "Helsingør"        
[19] "Herlev"            "Hillerød"          "Høje Taastrup"    
[22] "Hørsholm"          "Hvidovre"          "Ishøj"            
[25] "København"         "Lyngby-Taarbæk"    "Rødovre"          
[28] "Rudersdal"         "Tårnby"            "Vallensbæk"       
[31] "Århus"             "Favrskov"          "Hedensted"        
[34] "Herning"           "Holstebro"         "Horsens"          
[37] "Ikast-Brande"      "Lemvig"            "Norddjurs"        
[40] "Odder"             "Randers"           "Ringkøbing-Skjern"
[43] "Samsø"             "Silkeborg"         "Skanderborg"      
[46] "Skive"             "Struer"            "Syddjurs"         
[49] "Viborg"            "Aalborg"           "Brønderslev"      
[52] "Frederikshavn"     "Hjørring"          "Jammerbugt"       
[55] "Læsø"              "Mariagerfjord"     "Morsø"            
[58] "Rebild"            "Thisted"           "Vesthimmerland"   
[61] "Faxe"              "Greve"             "Guldborgsund"     
[64] "Holbæk"            "Kalundborg"        "Køge"             
[67] "Lejre"             "Lolland"           "Næstved"          
[70] "Odsherred"         "Ringsted"          "Roskilde"         
[73] "Slagelse"          "Solrød"            "Sorø"             
 [ reached getOption("max.print") -- omitted 24 entries ]

Instructions

The elections sf object should be already loaded in your environment.

• Load the cartogram package.
• Filter out the Danske Folkeparti votes from your elections dataset, creating a DF object
• Plot total electorate over municipality area for year 2015 in the DF data. Deviation from a straight line shows the degree of misrepresentation.
• Create a cartogram scaling to the pct_vote2015 column.
• Check that the DF voter population is proportional to the area.
• Plot the pct_vote2015 percentage on the cartogram. Notice how some areas have relatively shrunk or grown.

# load library
library(cartogram)

# Filter out Danske Folkeparti
DF <- elections %>%
    filter(grepl("^O", Party))
# Check the spread of votes and municipality area
plot(DF$pct_vote2015, st_area(DF, byid = TRUE), xlab = "Vote %", ylab = "Area (m2)",
    main = "Dansk Folkeparti fraction per municipality area")


# Make a cartogram, scaling the area to the percentage of SD voters
DF2015 <- cartogram_cont(DF, "pct_vote2015")

# Check the linearity of the SD voters percentage per municipality plot
plot(DF2015$pct_vote2015, st_area(DF2015, byid = TRUE))

Copacetic cartogram! Now try to rerun the cartogram for the Social Democrats in 2015 and create a visual for both parties’ turnout and total electorate in 2015.

library(cartogram)
# Let's look at Social Democrats in 2015
SD <- elections %>%
    filter(grepl("^A", Party))

# Make a cartogram, scaling the area to the total number of votes cast in 2015
SD2015 <- cartogram_cont(SD, "sum2015")

# Now check the linearity of the total voters per municipality cartogram as
# opposed to the reality
plot(SD$sum2015, st_area(SD, byid = TRUE))  # reality
plot(SD2015$sum2015, st_area(SD2015, byid = TRUE))  # cartogram


# Make a adjusted map of the 2015 SD and DF voters
plot(SD2015$geometry, col = "pink", main = "% of Social Democrat votes across DK in 2015")
plot(DF2015$geometry, col = "lightblue", main = "% of Danske Folkeparti votes across DK in 2015")

Instructions I - defining neighbors

• Load the elections spatial dataset with attributes
• Consider simplifying the boundaries if the data is too heavy for your computer and takes long to visualise
• Load the spdep library and create nb object of neighbors using queen adjacency
• Pass elections to poly2nb() to find the neighbors of each municipality polygon. Assign to nb.
• Get the center points of each municipality by passing elections to st_centroid and then to st_coordinates(). Assign to mun_centers.
• Update the basic map of the DK municipalities by adding the connections.
  • In the second plot call pass nb and mun_centers.
  • Also pass add = TRUE to add to the existing plot rather than starting a new one.

# Reload the data if needed elections
plot(elections$geometry)

# Consider simplifying (don't go too high)
mun_sm <- st_cast(st_simplify(mun, dTolerance = 250), to = "MULTIPOLYGON")
plot(mun_sm$geometry)
length(st_is_valid(mun_sm$geometry))

[1] 99

# Use the spdep package
library(spdep)

# Make neighbor list following queen adjacency
nb <- poly2nb(mun_sm$geometry)
nb

Neighbour list object:
Number of regions: 99 
Number of nonzero links: 344 
Percentage nonzero weights: 3.509846 
Average number of links: 3.474747 
11 regions with no links:
4 6 10 43 55 57 63 68 79 84 89
16 disjoint connected subgraphs

# Get center points of each municipality
mun_centers <- st_coordinates(st_centroid(mun_sm$geometry))

# Show the connections
plot(mun_sm$geometry)
plot(nb, mun_centers, col = "red", add = TRUE)

Instructions II - Moran’s I

Now that your neighbors are determined and centroids are computed, let’s continue with the Moran’s I statistic

• Create a subset with municipalities for O.Danske Folkeparti
• Feed the pct_2011 vector into moran.test().
  • moran.test() needs a weighted version of the nb object which you get by calling nb2listw().
  • After you specify your neighbor nb object you should define the weights style = "W". Here, style = "W" indicates that the weights for each spatial unit are standardized to sum to 1 (this is known as row standardization). For example, municipality 1 has 3 neighbors, and each of those neighbors will have weights of 1/3. This allows for comparability between areas with different numbers of neighbors.
  • You will need another argument in both spatial weights and at the level of the test. zero.policy= TRUE deals with situations when an area has no neighbors based on your definition of neighbor (many islands in Denmark). When this happens and you don’t include zero.policy= TRUE, you’ll get an error.
  • Run the test against the theoretical distribution of Moran’s I statistic. Find the p-value. Can you reject the null hypothesis of no spatial correlation?
• Inspect a map of pct_2011.
• Run another Moran I statistic test, this time on Social Democrats.
  • Use 999 Monte-Carlo iterations via moran.mc().
  • The first two arguments are the same as for moran.test().
  • You also need to pass the argument nsim = 999.
  • Note the p-value. Can you reject the null hypothesis this time?

# Let's look at Danske Folkeparti in 2015
DF <- elections %>% 
  _____(____)

# Run a Moran I test test on 2015 DF vote
moran.test(DF$_________, 
           nb2listw(____, style = "W",zero.policy=TRUE),
           zero.policy=TRUE)

# Run a Moran I test test on 2011 DF vote
moran.test(DF$________, 
           nb2listw(___, style = "W",zero.policy=TRUE),
           zero.policy=TRUE)

# Do a Monte Carlo simulation to get a more reliable p-value
moran.mc(DF$_________,
         ________(____, zero.policy=TRUE),
         zero.policy=TRUE, nsim = 999)

    Moran I test under randomisation

data:  DF$pct_vote2015  
weights: nb2listw(nb, style = "W", zero.policy = TRUE)  
n reduced by no-neighbour observations  

Moran I statistic standard deviate = 6.6385, p-value = 1.584e-11
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.532756868      -0.011494253       0.006721352 

    Moran I test under randomisation

data:  DF$pct_vote2011  
weights: nb2listw(nb, style = "W", zero.policy = TRUE)  
n reduced by no-neighbour observations  

Moran I statistic standard deviate = 6.2079, p-value = 2.685e-10
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.499144047      -0.011494253       0.006766158 

    Monte-Carlo simulation of Moran I

data:  DF$pct_vote2015 
weights: nb2listw(nb, zero.policy = TRUE)  
number of simulations + 1: 1000 

statistic = 0.53276, observed rank = 1000, p-value = 0.001
alternative hypothesis: greater

Marvelous Moran Testing! When I ran the examples, the p-value was around 1.584e-11 in 2015 and 2.685e-10 in 2011 Moran tests, showing significant spatial correlation. In Monte Carlo simulation, the p-value was around 0.001, so I did find some significant spatial (positive) correlation.

Repeat the same test for Social Democrats

# Let's look at Social Democrats

# Run a Moran I test on percentage of SD turnout in 2011

# Run a Moran I test on percent of SD turnout in 2015

# Do a Monte Carlo simulation to get a better p-value

# Do a Monte Carlo simulation to get a better p-value

    Moran I test under randomisation

data:  SD$pct_vote2011  
weights: nb2listw(nb, style = "W", zero.policy = TRUE)  
n reduced by no-neighbour observations  

Moran I statistic standard deviate = 5.6495, p-value = 8.044e-09
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.450511071      -0.011494253       0.006687554 

    Moran I test under randomisation

data:  SD$pct_vote2015  
weights: nb2listw(nb, style = "W", zero.policy = TRUE)  
n reduced by no-neighbour observations  

Moran I statistic standard deviate = 5.0148, p-value = 2.654e-07
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.40141853       -0.01149425        0.00677959 

    Monte-Carlo simulation of Moran I

data:  SD$pct_vote2011 
weights: nb2listw(nb, zero.policy = TRUE)  
number of simulations + 1: 1000 

statistic = 0.45051, observed rank = 1000, p-value = 0.001
alternative hypothesis: greater

    Monte-Carlo simulation of Moran I

data:  SD$pct_vote2015 
weights: nb2listw(nb, zero.policy = TRUE)  
number of simulations + 1: 1000 

statistic = 0.40142, observed rank = 1000, p-value = 0.001
alternative hypothesis: greater

Phenomenal political testing. Social Democrats also show positive correlation. p-value in Moran I test is was around 8.044e-09 in 2011 results and 2.654e-07 in 2015 results, and thus significant result. In Monte Carlo simulation, the p-value was around 0.24, suggesting there is significant (positive) spatial correlation.

Connect the nearest places (islands)

# Consider simplifying (don't go too high)
mun_sm <- st_cast(st_simplify(mun, dTolerance = 250), to = "MULTIPOLYGON")
plot(mun_sm$geometry)

# Use the spdep package
library(spdep)

# Get center points of each municipality
mun_centers <- st_centroid(mun_sm$geometry, of_largest_polygon = TRUE)

# Make neighbor list from neighbours at 100km distance
nb_100 <- dnearneigh(mun_centers, 0, 1e+05)
plot(mun_sm$geometry)
plot(nb_100, mun_centers, col = "red", add = TRUE)

# Make neighbor list from neighbours at 50km distance
nb_50 <- dnearneigh(mun_centers, 0, 50000)
plot(mun_sm$geometry)
plot(nb_50, mun_centers, col = "blue", add = TRUE)
title(main = "Neighbours within 50 km distance")